%% US Public Debt and Safe Asset Market Power
%% Jason Choi, Rishabh Kirpalani, and Diego Perez
%% Nov 24, 2024

%% Solve Monopoly Equilibrium

%----------------------------------------------------------------
% 0. Housekeeping
%----------------------------------------------------------------

close all

%----------------------------------------------------------------
% 1. Defining variables
%----------------------------------------------------------------

// Endogenous Variables
var b rb rkstar rk krw_star kus_star krw kus kstar k wstar w c_rw c_us drdb spread vrw vus dMrwdb by;

// Exogenous Variables
var nnu oomega A Astar eeta;

// Shocks
varexo eps_nnu eps_oomega eps_A eps_eeta;

// Parameters
parameters ggamma bbeta llambda aalpha Astarbar Abar iiota iiota_star ddelta_rw ddelta_us
  nnu_bar oomega_bar rrho_nnu rrho_oomega ssigma_nnu ssigma_oomega rrho_A ssigma_A rrho_Astar ssigma_Astar
  b_me rb_me rkstar_me rk_me krw_star_me kus_star_me krw_me kus_me kstar_me k_me
  capKstar_me capK_me wstar_me w_me crw_me cus_me spread_me vrw_me vus_me eeta_bar ssigma_eeta beta_alpha beta_beta beta_mean beta_var beta_sd;

%----------------------------------------------------------------
% 2. Calibration
%----------------------------------------------------------------

% // Parameters
disp('% Epsilon = 2.2, Lambda = 1') 
ggamma = 2;
bbeta = 0.9886;
% eeta = 0.545;
llambda = 1;
aalpha = 0.3;
Astarbar = 0.9254;
Abar = 0.8154;
iiota = 0.9070;
iiota_star = 0.7939;
ddelta_rw = 0.1;
ddelta_us = 0.1;
nnu_bar = 0.0048;
oomega_bar = 0.0025;
rrho_nnu = 0.99;
ssigma_nnu = 0.01;
rrho_oomega = 0.95;
ssigma_oomega = 0.48;
rrho_A = 0.95;
ssigma_A = 0.02;
rrho_Astar = rrho_A;
ssigma_Astar = ssigma_A;

eeta_bar = 0.545;

beta_alpha = 0.001;
beta_beta = beta_alpha/eeta_bar - beta_alpha;
beta_mean = beta_alpha/(beta_alpha+beta_beta);
beta_var = (beta_alpha*beta_beta)/((beta_alpha+beta_beta)^2*(beta_alpha+beta_beta+1));
beta_sd = sqrt(beta_var);
ssigma_eeta = beta_sd;

% Analytic Steady State (Monopoly Equilibrium)
b_me = (oomega_bar/(nnu_bar*eeta_bar))^(1/(eeta_bar-1-llambda));
rb_me = 1/bbeta - nnu_bar*(b_me)^(eeta_bar-1) - 1;
rkstar_me = 1/bbeta + ddelta_rw - 1;
rk_me = 1/bbeta + ddelta_us - 1;
krw_star_me = ((aalpha*(1-iiota_star)*Astarbar*((1/bbeta+ddelta_rw-1)/(aalpha*iiota_star*Astarbar))^((aalpha*(1-iiota_star)-1)/(aalpha*(1-iiota_star))))/(1/bbeta+ddelta_us-1))^((aalpha*(1-iiota_star))/(1-aalpha));
kus_star_me = ((1/bbeta+ddelta_rw-1)/(aalpha*iiota_star*Astarbar))^(1/(aalpha*(1-iiota_star)))*krw_star_me^((1-iiota_star*aalpha)/(aalpha*(1-iiota_star)));
krw_me = ((aalpha*iiota*Abar*((1/bbeta+ddelta_us-1)/(aalpha*(1-iiota)*Abar))^((aalpha*iiota-1)/(aalpha*iiota)))/(1/bbeta+ddelta_rw-1))^((aalpha*iiota)/(1-aalpha));
kus_me = ((1/bbeta+ddelta_us-1)/(aalpha*(1-iiota)*Abar))^(1/(aalpha*iiota))*krw_me^((1-(1-iiota)*aalpha)/(aalpha*iiota));
kstar_me = krw_star_me + krw_me;
k_me = kus_star_me + kus_me;
capKstar_me = krw_star_me^iiota_star*kus_star_me^(1-iiota_star);
capK_me = krw_me^(1-iiota)*kus_me^iiota;
wstar_me = Astarbar*(1-aalpha)*(capKstar_me)^aalpha;
w_me = Abar*(1-aalpha)*(capK_me)^aalpha;
crw_me = wstar_me + (rkstar_me-ddelta_rw)*kstar_me + nnu_bar/eeta_bar*(b_me)^eeta_bar + rb_me*b_me;
cus_me = w_me + (rk_me-ddelta_us)*k_me - oomega_bar/(1+llambda)*(b_me)^(1+llambda) - rb_me*(b_me);
drdb_me = -nnu_bar*(eeta_bar-1)*b_me^(eeta_bar-2);
dMrwdb_me = 0;
nnu_me = nnu_bar;
oomega_me = oomega_bar;
A_me = Abar;
Astar_me = Astarbar;
spread_me = (rk_me-ddelta_us-rb_me);
vrw_me = crw_me^(1-ggamma)/(1-ggamma)/(1-bbeta);
vus_me = cus_me^(1-ggamma)/(1-ggamma)/(1-bbeta);
by_me = b_me/(A_me*(kus_me^iiota*krw_me^(1-iiota))^aalpha);
eeta_me = eeta_bar;

%----------------------------------------------------------------
% 3. Model
%----------------------------------------------------------------

model;

c_rw^(-ggamma) = bbeta*c_rw(+1)^(-ggamma)*(nnu*b^(eeta-1)+1+rb);
c_rw^(-ggamma) = bbeta*c_rw(+1)^(-ggamma)*(1-ddelta_rw+rkstar);
c_rw + kstar + b = wstar + (1-ddelta_rw+rkstar(-1))*kstar(-1) + nnu(-1)/eeta(-1)*(b(-1))^eeta(-1) + (1+rb(-1))*b(-1);

c_us^(-ggamma) = bbeta*c_us(+1)^(-ggamma)*(oomega*b^llambda+1+rb+drdb*b);
c_us^(-ggamma) = bbeta*c_us(+1)^(-ggamma)*(1-ddelta_us+rk);
c_us + k - b = w + (1-ddelta_us+rk(-1))*k(-1) - oomega(-1)/(1+llambda)*(b(-1))^(1+llambda) - (1+rb(-1))*b(-1);

rk = Astar*aalpha*(1-iiota_star)*krw_star^(aalpha*iiota_star)*kus_star^(aalpha*(1-iiota_star)-1);
rkstar = Astar*aalpha*iiota_star*krw_star^(aalpha*iiota_star-1)*kus_star^(aalpha*(1-iiota_star));
rk = A*aalpha*iiota*krw^(aalpha*(1-iiota))*kus^(aalpha*iiota-1);
rkstar = A*aalpha*(1-iiota)*krw^(aalpha*(1-iiota)-1)*kus^(aalpha*iiota);
wstar = Astar*(1-aalpha)*krw_star^(aalpha*iiota_star)*kus_star^(aalpha*(1-iiota_star));
w = A*(1-aalpha)*krw^(aalpha*(1-iiota))*kus^(aalpha*iiota);

0 = -(drdb+nnu*(eeta-1)*b^(eeta-2))*(c_rw(+1)/c_rw)^(-ggamma)+(1+rb+nnu*b^(eeta-1))*dMrwdb;
0 = dMrwdb*(rkstar+1-ddelta_rw);

k = kus + kus_star;
kstar = krw + krw_star;

eeta = eeta_bar + ssigma_eeta*eps_eeta;

log(nnu) = (1-rrho_nnu)*log(nnu_bar) + rrho_nnu*log(nnu(-1)) + ssigma_nnu*eps_nnu;
log(oomega) = (1-rrho_oomega)*log(oomega_bar) + rrho_oomega*log(oomega(-1)) + ssigma_oomega*eps_oomega;
log(A) = (1-rrho_A)*log(Abar) + rrho_A*log(A(-1)) + ssigma_A*eps_A;
log(Astar) = (1-rrho_Astar)*log(Astarbar) + rrho_Astar*log(Astar(-1)) + ssigma_Astar*eps_A;

spread = (rk-ddelta_us-rb);

vrw = c_rw^(1-ggamma)/(1-ggamma) + bbeta*vrw(+1);
vus = c_us^(1-ggamma)/(1-ggamma) + bbeta*vus(+1);

by = b/(A*(kus^iiota*krw^(1-iiota))^aalpha);

end;

%----------------------------------------------------------------
% 4. Computation
%----------------------------------------------------------------

initval;
  b = b_me;
  rb = rb_me;
  rkstar = rkstar_me;
  rk = rk_me;
  krw_star = krw_star_me;
  kus_star = kus_star_me;
  krw = krw_me;
  kus = kus_me;
  kstar = kstar_me;
  k = k_me;
  wstar = wstar_me;
  w = w_me;
  c_rw = crw_me;
  c_us = cus_me;
  drdb = drdb_me;
  dMrwdb = dMrwdb_me;
  nnu = nnu_me;
  oomega = oomega_me;
  spread = spread_me;
  A = A_me;
  Astar = Astar_me;
  vrw = vrw_me;
  vus = vus_me;
  by = by_me;
  eeta = eeta_me;
end;

resid;
check;

shocks;
    var eps_nnu = 1;
    var eps_A = 1;
    var eps_oomega = 1;
    var eps_eeta = 1;
end;

set_dynare_seed('default');
stoch_simul(order=2,nograph,noprint,periods=10000,pruning);

%----------------------------------------------------------------
% 5. Generate moments
%----------------------------------------------------------------

spread_path = (rk-ddelta_us-rb)*100;
var_b = var(b);
var_sp = var(spread_path);
auto_b = autocorr(b);
auto_sp = autocorr(spread_path);
cost = oomega./(1+llambda).*(b).^(1+llambda);
benefit = (eeta-1).*b.*(log(b)-nnu);
cost_p = oomega.*(b).^(llambda);
benefit_p = (eeta-1)./b;
corr_pq = corr(spread_path,b);
var_by = var(by);
auto_by = autocorr(by);
corr_pq_by = corr(spread_path,by);

moments = [mean(by) mean(spread_path) var_by var_sp corr_pq_by auto_by(2) auto_sp(2)]';
data_mom = [0.41 0.62 0.03 0.086 -0.56 0.96 0.70]';
rowNames = {'Mean b/y','Mean sp','Var b/y','Var sp','Corr (b/y,sp)','Autocorr b/y','Autocorr sp'};
colNames = {'Model Moments','Data Moments'};
TableA0 = array2table([moments data_mom],'RowNames',rowNames,'VariableNames',colNames)

deficit_by = cost(1:end-1) + by(2:end) - (1+rb(1:end-1)).*by(1:end-1);
ca_by = -(by(2:end)-by(1:end-1)) + kus_star(2:end)-kus_star(1:end-1) - (krw(2:end)-krw(1:end-1));
nfa_by = - by + kus_star - krw;
var_ca_by = var(ca_by);
var_nfa_by = var(nfa_by);
var_deficit_by = var(deficit_by);
corr_nfa_by = corr(nfa_by,by);
corr_ca_def_by = corr(ca_by,deficit_by);

moments = [mean(rb) var(rb)*100 var_ca_by*100 var_nfa_by var_deficit_by corr_nfa_by corr_ca_def_by]';
data_mom = [0.0053 0.00097*100 0.00033*100 0.035 0.002 -0.654 -0.207]';
rowNames = {'Mean rb (Targeted)','Var rb(x100)','Var CA(x100)','Var NFA','Var Deficit','Corr(NFA,b)','Corr(CA,deficit)'};
colNames = {'Model Moments','Data Moments'};
TableA0 = array2table([moments data_mom],'RowNames',rowNames,'VariableNames',colNames)

%----------------------------------------------------------------
% 6. Save moments for welfare calculations
%----------------------------------------------------------------
b_pos = strmatch('b',M_.endo_names,'exact');
rb_pos = strmatch('rb',M_.endo_names,'exact');
rkstar_pos = strmatch('rkstar',M_.endo_names,'exact');
rk_pos = strmatch('rk',M_.endo_names,'exact');
krw_star_pos = strmatch('krw_star',M_.endo_names,'exact');
kus_star_pos = strmatch('kus_star',M_.endo_names,'exact');
krw_pos = strmatch('krw',M_.endo_names,'exact');
kus_pos = strmatch('kus',M_.endo_names,'exact');
kstar_pos = strmatch('kstar',M_.endo_names,'exact');
k_pos = strmatch('k',M_.endo_names,'exact');
wstar_pos = strmatch('wstar',M_.endo_names,'exact');
w_pos = strmatch('w',M_.endo_names,'exact');
crw_pos = strmatch('c_rw',M_.endo_names,'exact');
cus_pos = strmatch('c_us',M_.endo_names,'exact');
drdb_pos = strmatch('drdb',M_.endo_names,'exact');
dMrwdb_pos = strmatch('dMrwdb',M_.endo_names,'exact');
vrw_pos = strmatch('vrw',M_.endo_names,'exact');
vus_pos = strmatch('vus',M_.endo_names,'exact');
nnu_pos = strmatch('nnu',M_.endo_names,'exact');
oomega_pos = strmatch('oomega',M_.endo_names,'exact');
A_pos = strmatch('A',M_.endo_names,'exact');
Astar_pos = strmatch('Astar',M_.endo_names,'exact');
z_pos = strmatch('z',M_.endo_names,'exact');
by_pos = strmatch('by',M_.endo_names,'exact');
eeta_pos = strmatch('eeta',M_.endo_names,'exact');

shock_matrix = randn(size(oo_.endo_simul,2),M_.exo_nbr); %create shock matrix with number of time periods in rows

shock_matrix(:,4) = (betarnd(beta_alpha,beta_beta,size(oo_.endo_simul,2),1) - eeta_bar)/beta_sd; %create shock matrix with number of time periods in rows

sim_y_forward = simult_(M_,options_,oo_.mean,oo_.dr,shock_matrix,options_.order);

z = sim_y_forward(oomega_pos,:); 
b_path = sim_y_forward(by_pos,:);
FVD = sim_y_forward(eeta_pos,:)<eeta_bar;
eeta_path = sim_y_forward(eeta_pos,:);
spread_path = (sim_y_forward(rk_pos,:)-ddelta_us-sim_y_forward(rb_pos,:))*100;

spread_path(b_path<0)=[];
z(b_path<0)=[];
FVD(b_path<0)=[];
eeta_path(b_path<0)=[];
b_path(b_path<0)=[];

b_path(spread_path<0)=[];
z(spread_path<0)=[];
FVD(spread_path<0)=[];
eeta_path(spread_path<0)=[];
spread_path(spread_path<0)=[];

XX = log(b_path);
YY = spread_path;

XXb = [log(b_path)'];
XX = [log(b_path)' FVD'];
YY = spread_path;

XXp = [log(b_path)' FVD' log(b_path)'.*FVD'];

table = fitlm(XXp,YY)
low_beta = table2array(table.Coefficients(2,1));
high_beta = table2array(table.Coefficients(2,1))+table2array(table.Coefficients(4,1));
low_elast = 1/(low_beta/mean(spread_path));
high_elast = 1/(high_beta/mean(spread_path));

model_mom = [high_elast' low_elast']';
data_mom = [-1.45 -3.23]';
rowNames = {'High Vol Elasticity','Low Vol Elasticity'};
colNames = {'Model Moments','Data Moments'};
TableA0 = array2table([model_mom data_mom],'RowNames',rowNames,'VariableNames',colNames)

[coef, stderr, pval] = tsls(YY(2:end)',[ones(length(b_path)-1,1) log(b_path(2:end))' FVD(2:end)' log(b_path(2:end))'.*FVD(2:end)'],[ones(length(b_path)-1,1) FVD(2:end)' z(2:end)' log(b_path(1:end-1))' z(2:end)'.*FVD(2:end)' log(b_path(1:end-1))'.*FVD(1:end-1)' ])

low_beta = coef(2);
high_beta = coef(2)+coef(4);
low_elast = 1/(low_beta/mean(spread_path));
high_elast = 1/(high_beta/mean(spread_path));

model_mom = [high_elast' low_elast']';
data_mom = [-1.45 -3.23]';
rowNames = {'High Vol Elasticity','Low Vol Elasticity'};
colNames = {'Model Moments','Data Moments'};
TableA0 = array2table([model_mom data_mom],'RowNames',rowNames,'VariableNames',colNames)

writematrix([XXp YY' z' eeta_path'],'simulated_me.xlsx')

%----------------------------------------------------------------
% 6. Calculate welfare from transition
%---------------------------------------------------------------

b_me_eta = mean(b);
rb_me_eta = mean(rb);
spread_me_eta = mean(spread);

oo_me_eta = oo_;
M_me_eta = M_;
options_me_eta = options_;

save me_eta_save oo_me_eta vus_me vrw_me cus_me crw_me M_me_eta options_me_eta b_me_eta spread_me_eta rb_me_eta;

fileID = fopen('beta_reg.tex','w');
ddata = [table2array(table.Coefficients(2,1)) coef(2) table2array(table.Coefficients(2,2)) stderr(2) table2array(table.Coefficients(3,1)) coef(3) table2array(table.Coefficients(3,2)) stderr(3) table2array(table.Coefficients(4,1)) coef(3) table2array(table.Coefficients(4,2)) stderr(3)]';
fprintf(fileID,'\\begin{tabular}{l>{\\centering}p{1.5cm}>{\\centering}p{1.5cm}>{\\centering}p{1.5cm}>{\\centering}p{1.5cm}} \n \\toprule \n \\multicolumn{1}{c}{} & \\multicolumn{2}{c}{Data} & \\multicolumn{2}{c}{Simulations}\\tabularnewline \n VARIABLES & OLS & IV & OLS & IV\\tabularnewline \n \\midrule \n Log(debt/gdp) & -0.23{*}{*}{*} \\vspace{0.5em} \n  & -0.16{*}{*}{*} & %3.2f{*}{*}{*} & %3.2f{*}{*}{*}\\tabularnewline \n  & (0.04) \\vspace{0.5em} \n  & (0.05) & (%3.2f) & (%3.2f)\\tabularnewline \n Demand Rotator $z$ & -0.31{*}{*}{*} \\vspace{0.5em} \n  & -0.12 & %3.2f{*}{*}{*} & %3.2f{*}{*}{*}\\tabularnewline \n  & (0.07) \\vspace{0.5em} \n  & (0.09) & (%3.2f) & (%3.2f)\\tabularnewline \n $z$ $\\times$ Log(debt/gdp) & -0.34{*}{*}{*} \\vspace{0.5em} \n  & -0.18{*}{*} & %3.2f{*}{*}{*} & -%3.2f{*}{*}{*}\\tabularnewline \n  & (0.06) \\vspace{0.5em} \n  & (0.08) & (%3.2f) & (%3.2f)\\tabularnewline \n US VIX & 0.02{*}{*}{*} \\vspace{0.5em} \n  & 0.02{*}{*}{*} &  & \\tabularnewline \n  & (0.00) \\vspace{0.5em} \n  & (0.00) &  & \\tabularnewline \n Yield curve slope & -0.06{*}{*}{*} \\vspace{0.5em} \n  & -0.07{*}{*}{*} &  & \\tabularnewline \n  & (0.01) \\vspace{0.5em} \n  & (0.01) &  & \\tabularnewline \n Constant & 0.11{*}{*} \\vspace{0.5em} \n  & 0.17{*}{*}{*} & 0.20{*}{*}{*} & 0.21{*}{*}{*}\\tabularnewline \n  & (0.06) \\vspace{0.5em} \n  & (0.07) & (0.01) & (0.01)\\tabularnewline \n \\midrule \n Observations & 339 \\vspace{0.5em} \n & 339 & 100,000 & 100,000\\tabularnewline \n R-squared & 0.53 \\vspace{0.5em} \n &  & 0.63 & \\tabularnewline \n Demand elasticity, high foreign vol & 1.09 \\vspace{0.5em} \n & 1.82 & 2.22 & 1.54\\tabularnewline \n Demand elasticity, low foreign vol & 2.68 \\vspace{0.5em} \n & 3.78 & 3.25 & 3.62\\tabularnewline \n Markup, high foreign vol & 0.92 \\vspace{0.5em} \n & 0.55 & 0.45 & 0.65\\tabularnewline \n Markup, low foreign vol & 0.37 \\vspace{0.5em} \n & 0.26 & 0.31 & 0.28\\tabularnewline \n \\bottomrule \n \\end{tabular}',ddata);
fclose(fileID);

fileID = fopen('beta_conduct.tex','w');
ddata = [table2array(table.Coefficients(2,1)) coef(2) table2array(table.Coefficients(2,2)) stderr(2) table2array(table.Coefficients(3,1)) coef(3) table2array(table.Coefficients(3,2)) stderr(3) table2array(table.Coefficients(4,1)) coef(3) table2array(table.Coefficients(4,2)) stderr(3)]';
fprintf(fileID,'\\begin{tabular}{>{\\raggedright}p{2cm}>{\\centering}p{1.3cm}>{\\centering}p{1.3cm}>{\\centering}p{1.3cm}>{\\centering}p{1.3cm}>{\\centering}p{1.3cm}>{\\centering}p{1.3cm}>{\\centering}p{1.3cm}>{\\centering}p{1.3cm}>{\\centering}p{1.3cm}} \n \\toprule \n Cost elasticity & {\\small{}$\\lambda=0$} & {\\small{}$\\lambda=0.5$} & {\\small{}$\\lambda=1$} & {\\small{}$\\lambda=2$} & {\\small{}$\\lambda=3$} & {\\small{}$\\lambda=4$} & {\\small{}$\\lambda=10$} & {\\small{}$\\lambda=20$} & {\\small{}$\\lambda=30$}\\tabularnewline \n \\midrule \n \\emph{a. Monopoly} \\vspace{0.5em} \n &  &  &  &  &  &  &  &  & \\vspace{0.5em} \n \\tabularnewline \n & -15.73 & -12.65 & -10.14 & -7.28 & -5.84 & -5.00 & -3.25 & -2.60 & -2.38 \\vspace{0.5em} \n \\tabularnewline \n & (0.0) & (0.0) & (0.0) & (0.0) & (0.0) & (0.0) & (0.0) & (0.0) & (0.0) \\vspace{0.5em} \n \\tabularnewline \n \\emph{b. Price-taking} \\vspace{0.5em} \n &  &  &  &  &  &  &  &  & \\vspace{0.5em} \n \\tabularnewline \n & 5.11 & 10.03 & 12.93 & 15.57 & 16.60 & 17.10 & 17.82 & 17.98 & 18.02 \\vspace{0.5em} \n \\tabularnewline \n & (0.0) & (0.0) & (0.0) & (0.0) & (0.0) & (0.0) & (0.0) & (0.0) & (0.0) \\vspace{0.5em} \n \\tabularnewline \n \\bottomrule \n \\end{tabular}');
fclose(fileID);